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Abstract 

I report results of a density matrix renormalization group (DMRG) study 

of a model for the two dimensional spin-gapped system CaV^g. This study 

represents the first time that DMRG has been used to study a two dimensional 

system on large lattices, in this case as large as 24 x 11, allowing extrapolation 

to the thermodynamic limit. I present a substantial improvement to the 

DMRG algorithms which makes these calculations feasible. 
PACS Numbers: 75.10.-b., 75.10.Jm, 75.40. Gb 
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Since the discovery of the high temperature superconductivity, condensed matter physi- 
cists have searched for other two dimensional systems with exotic "spin liquid" ground 
states. Thus considerable excitement accompanied the recent discovery of CaX^Og, a two 
dimensional, frustrated S = 1/2 Heisenberg spin system, with a substantial spin gap [|TJ. 
This system has been modeled by a depleted square lattice Heisenberg model, with both 
nearest and next-nearest exchange interactions ||. It consists of a square lattice with 1/5 
of the spins missing, as shown in Figure 1, and has been called a CAVO lattice ||. It is 
believed that the superexchange is mediated by out-of-plane oxygen atoms, resulting in a 
very large next-nearest exchange: J' ~ J/2 |Q. This frustrating interaction helps stabilize 
the spin-gapped state against Neel order. 

One can think of the ground state of this system as a "plaquette resonating valence 
bond" state 0: no phase transition is expected if the interactions between plaquettes are 
adiabatically removed, and the ground state of a single plaquette is perfectly described by a 
resonating valence bond (RVB) variational ansatz (in the absence of frustration). In addi- 
tion, in the weakly- interacting plaquette limit, pairs of holes bind on plaquettes, suggesting 
the possibility of a superconducting ground state upon doping. The system is reminiscent of 
ladder systems with even numbers of legs. In that case, somewhat more complicated RVB 
states have been useful in describing the spin liquid ground states of undoped ladders ||||, 
and upon doping, strong pairing correlations are observed numerically |J. 

A number of theoretical and numerical treatments have been performed on this model in 
the last year B^HTTI]. Troyer, Kontani and Ueda made the most reliable determination of the 



unfrustrated phase diagram using a quantum Monte Carlo loop algorithm fL0| . An important 
conclusion of this study was that simple 1/5 depletion of the isotropic square lattice, without 
frustration, does not destroy Neel order. A spin liquid ground state was found when the 
couplings within a plaquette were about 10% greater than between plaquettes. This result 
contradicted earlier (non-loop) quantum Monte Carlo calculations on smaller systems 
It was not possible to include frustration because of sign problems. Gelfand, et. al. |3j 
applied series expansion techniques and were able to study the frustrated and unfrustrated 



systems. Their results were in agreement with Troyer, et. al. for the unfrustrated system. 
They concluded that a next-nearest neighbor interaction J' = J/2 was consistent with 
experimental results. 

I present results here from density matrix renormalization group (DMRG) calculations 
|T^| for the spin gap of the frustrated CAVO lattice. The results are in agreement with those 



of Gelfand, et. al., and, in fact, when extrapolated to the thermodynamic limit, appear to 
be more accurate. Although DMRG is usually much more accurate than other numerical 
techniques for large one dimensional systems, these are the first reliable results for systems 
wide enough to be considered two dimensional — up to 24 x 11. These calculations are feasible 
because of an important improvement to the DMRG algorithms, which I present here, which 
increases the speed of the calculations by up to two orders of magnitude. 

The improvement to DMRG involves keeping track of the wavefunction from step to 
step. The step referred to here is the process of adding a site to a block and requires the di- 



agonalization of a superblock configuration of two blocks and two sites [|T2|| . In each DMRG 
step, an iterative sparse matrix algorithm, such as the Davidson method, is used to find 
the ground state of the superblock. In the original formulation of DMRG, no starting point 
for the Davidson procedure was specified. To ensure that the DMRG procedure is always 
stable and convergent, the superblock ground state usually has to be determined to rather 
high accuracy. (One diagonalization which converges to a low-lying eigenstate other than 
the ground state ruins the accuracy of the entire DMRG calculation.) Consequently, a sub- 
stantial number of Davidson steps are necessary to converge to sufficient accuracy, typically 
40-100. The total calculation time is proportional to the average number of Davidson steps. 

If a very good initial guess is available for the Davidson procedure, the number of David- 
son steps can be reduced substantially. An ideal initial guess, for the case of the finite 
system DMRG algorithm, is the final wavefunction from the previous DMRG step. This 
wavefunction, however, is in a different basis, corresponding to a different superblock, but 
it can be transformed into the basis corresponding to the current superblock, as I describe 
below. Use of this transformation to obtain the initial state in a Davidson diagonalization 



can reduce the number of Davidson steps by one half, typically, assuming that one iterates 
Davidson until it converges to high accuracy. Use of this initial guess has an even more 
important advantage: it is not necessary to converge to high accuracy, since there is no 
danger of converging to an incorrect low-lying eigenstate. The initial guess not only has low 
energy, it approximately describes the correct eigenstate, as obtained in the previous step. 
In fact, the algorithm can be made completely stable even if the number of Davidson steps 
is restricted to two or three! Thus one saves a factor of 20-50 in the time required by the 
Davidson procedure. The overall speedup is somewhat reduced from this factor because the 
calculation time to perform other parts of the DMRG procedure, such as diagonalizing the 
density matrix, becomes significant. 

A DMRG step adds a site onto a block, constructing an appropriate basis for the new 
block. Let \ai) be the states of left block /, where I is the rightmost site of the block. Let 
\si) be the states of site I. Then the basis states for the new left block are given by 

|a I+ i) = L l+1 [s l+ i] ai+uai \ai) <g> |sj+i). (1) 



This notation is similar to that of Ostlund and Rommer jTJ]. The transformation matrix 
L l+1 [si + i] ai+lt0ll is a slightly rewritten form of the truncated matrix of density matrix eigen- 
vectors. The states of the right block \/3i+3) were formed at an earlier DMRG step in a 
similar fashion 

|A +3 >= E R l+3 [si + s]p l+3A+4 \s l+3 )^\(3 l+i ). (2) 

I do not assume any reflection symmetry for the lattice: the L and R matrices are indepen- 
dent. 

A superblock basis state is written in the form 

\aiSi +1 si +2 Pi +3 ) = K) ® \si+i) <S> \si +2 ) <8> | A+3>- (3) 
A superblock wavefunction is written in this basis as 

\^) = E i>(®iSi + iSi +2 (3i + 3)\atiSi + iSi + 2l3i +3 ). (4) 

aiSi +1 s l+2 (3i +3 



One needs to transform this wavefunction into the basis appropriate for the next DMRG 
step, \ai + iSi + 2Si +3 j3i +4 ) . The transformation is not exact, since there is a truncation in going 
from \aisi + i) to |cfy + i). However, the states \eti+i) are formed using the density matrix to 
be ideally adapted for representing so for the transformation of the wavefunction only, 
one can approximate 

\a l+1 ){a l+1 \ « 1. (5) 
With this approximation one readily obtains 

ll>(a l+1 Si +2 S l+3 P l+4 ) « Yl ^ +1 [ s i+l]a ;+ i,aXa^i + iSi +2 A+3)^ +3 [Si + 3]A + 3,/3 i+ 4- ( 6 ) 

O-lSl+lPl+Z 

The most efficient way to implement this transformation numerically is to first form the 
intermediate wavefunction 

i/j(a l+1 s l+2 (3i +3 ) = J2 ^ +1 [ s i+i]« ; +i,a ! ^(aiSi+is«+2A+3), (7) 
and then form the final result 

^(ai+lSj + 2Si+3/3i+4) = S ^( a l+l S l+2pl+3)R l+3 [si+3]p l+3 ,f3 l+ 4- ( 8 ) 

A+3 

In this form, the transformation requires very computer little time compared to other parts 
of the calculation. 

This transformation is used for one half of the DMRG steps, when a site is being added 
to the left block. An analogous transformation is used for adding a site to the right block. 

Implementing this transformation requires saving all the matrices L and R, which is 
ordinarily not done. The storage for these matrices is typically 20-30% of the storage required 
for the blocks themselves, so the extra storage is not a major concern. In an efficient DMRG 
implementation for a typical machine, such as a Cray or a workstation, both the blocks 
and the transformation matrices should be stored on disk. The calculations described here 
sometimes required more than a gigabyte of scratch disk storage, but never more than 80-90 
megabytes of RAM. 
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In many one dimensional systems, DMRG converges in one or two sweeps through the 
system. In quasi-one or two dimensional systems, the number of sweeps needed can easily 
grow to five to 10. Another important improvement in efficiency comes from gradually 
increasing the number of states kept per block as one performs the sweeps |14|]. In this 



case the calculation time is dominated by the last sweep, during which the number of 
Davidson steps per DMRG step can be constrained to be only two or three. Compared to 
a DMRG calculation keeping a constant number of states and without the wavefunction 
transformation, the speedup can be over two orders of magnitude. 

If the lattice is reflection-symmetric, an additional factor of two can be saved in both 
calculation time and memory, but a different wavefunction transformation is needed for the 
DMRG step where the superblock is symmetric. Also, it is possible to adopt somewhat 
similar methods to obtain good initial guesses for the wavefunction in the infinite system 
method. These techniques will be reported elsewhere |15| . 



The beginning of any DMRG calculation of a 2D system is a mapping of the 2D lattice 
onto a ID chain — basically, one must choose an order to traverse the sites. It is standard to 
use the scanline mapping — fix x, step through all values of y, then increment x, etc. — which 
has the advantage of keeping the blocks as contiguous as possible. In treating the CAVO 
lattice, it is more advantageous to modify this slightly so that all the sites in a plaquette are 
traversed in succession. This incorporates the fact that for the parameters I consider here, 
correlations within a plaquette are strongest. 

I consider the Heisenberg Hamiltonian 

H = J l3 £ S, • Sj (9) 

(id) 

defined on an L x x L y CAVO lattice with S — |. As shown in Fig. 1, I take all nearest- 
neighbor Jij to be identical, with value J\ — 1, setting the energy scale. All next-nearest- 
neighbor Jij are also identical, with value J2. All other Jy are zero. 

I have studied open systems, although putting periodic boundary conditions in the short 
direction y is not particularly difficult. Having open boundary conditions allows a variety 
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of sizes to be studied. We consider a set of systems of length 24 and width up to 11. Since 
interactions are strongest within plaquettes, only full plaquettes are included. We keep up 
to m — 600 states per block, with truncation error at worst about 10~ 5 . Typical errors in 
the total energy, for the larger systems, were less than about 10~ 3 . For systems of width 8 
and larger, an extrapolation to m — > oo was used, assuming an exponential fall off in the 
error in the energy as a function of m [|I(^,|l7[]. Corrections were at most about 10 -3 . For each 
system we calculate the ground state energies with quantum numbers S z = and S z = 1. 
The spin gap A is the difference in energies. The largest system, 24 x 11, took about 15 - 
20 hours of workstation time (rated at 135 SPECfp92) for one value of S z . 

Figure 2 shows some of the results, for various widths of the system L y . From the width 
7 data, we see that the spin gap is peaked at J 2 = 0.5, which also happens to be appropriate 
for CaX^Og. For larger values of J 2 it falls rapidly towards zero. (The gap for J 2 = 0.8 was 
consistent with zero, within uncharacteristically large error bars of about 0.05.) 

Finite size extrapolation is crucial to determine the spin gap for smaller values of J 2 . 
Excellent extrapolations can be obtained if one assumes the low lying spin excitations obey 
a relativistic dispersion relation 

A(k) 2 = A 2 + v 2 k 2 , (10) 

where A is the bulk gap and the velocity v corresponds to the speed of light. Lorentz 
invariant low energy excitations are common (but not universal) in gapped, one dimensional 
spin systems, reflecting Lorentz invariance of the corresponding nonlinear sigma model. 
Considerations of simple particle-in-a-box systems with various boundary conditions indicate 
that a generic boundary condition at the edge of an open system fixes the logarithmic 
derivative of the wavefunction ijj, j^/ip = const . This can be shown to imply that the 
lowest value of k allowed in a 1-D box of size L is given by 7r/ (L — a), where a depends on 
boundary effects but is independent of L. This leads to the following form for the g ap as a 
function of system size 

x9 .9 ^V 2 H 2 V 2 , , 
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I have found that for J 2 < ~ 0.3, we can set a = 0, and still obtain excellent fits. This 
is indicated by linear behavior when A 2 is plotted versus 1/Ly, with L x fixed. Results are 
shown in Figure 3 for typical values of J 2 . For J 2 = 0.5, the fits were poor for a = 0, and 
mediocre with a nonzero, because the data was slightly irregular. Gelfand, et.al. || observed 
that the gap minimum can move away from (ir, tt) for larger values of J 2 . This data supports 
that proposition for J 2 = 0.5. In the case of an incommensurate gap minimum, we would 
expect irregular behavior of the g function of L, as the value of k allowed by the 

lattice which is closest to the minimum would jump about as L increased. 

Using the fits shown, I corrected for the finite value of L x and obtained results for A in 
the thermodynamic limit. Figure 4 shows the results as a function of J r 2 . The extrapolation 
using Eq. (11) yielded imaginary gaps for J 2 = and 0.05, which we interpret to mean 
A = 0. The transition to a spin gapped state appears at J2 = 0.06(1). This result is 
in agreement with previous work indicating that the system with J2 = is close to the 
disordered phase. The value at J 2 = 0.5, A = 0.515(15), is somewhat lower than the series 
results of Gelfand, et.al., A = 0.57(3). However, the results are completely consistent if the 
shift in the gap minimum results in an overestimate in the series results by about 0.05, as 
Gelfand, et.al. suggest. 

I have presented the first results using DMRG on a two dimensional system for lattices 
wide enough to allow extrapolation to the thermodynamic limit. DMRG still is primarily 
a one dimensional technique, in that the accuracy falls off rapidly as the system's width 
increases. The CAVO system studied here was less difficult than many other two dimensional 
systems, both because of the existence of a gap and the depleted character of the lattice. 
Nevertheless, as DMRG and computational resources improve, I expect it will become a 
standard numerical technique for two dimensional systems, including doped fermion systems. 

I thank Naokazu Shibata and Rajiv Singh for helpful conversations. I acknowledge 
support from the Office of Naval Research under grant No. N00014-91-J-1143, and from the 
NSF under Grant No. DMR-9509945. Some of the calculations were performed at the San 
Diego Supercomputer Center. 
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FIGURES 

FIG. 1. A 12 x 7 open CAVO lattice. The solid lines are nearest-neighbor bonds with ex- 
change J\ , and the dotted lines are next-nearest-neighbor bonds with exchange Ji . Only complete 
plaquettes which fit within a 12 x 7 rectangle are retained. 



FIG. 2. Spin gap as a function of J2 for various widths L y , with L x = 24. 

FIG. 3. Gap squared as a function of L~ 2 . The solid lines are linear fits, excluding the L y = 4 
point for J' = 0.2. The dashed line is an alternative fit which includes an L y 3 term. 

FIG. 4. Spin gap extrapolated to the thermodynamic limit, as a function fo J2. Where not 
shown, errors are comparable to the point size. 
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